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We present spatially localized nonrotating and rotating (azimuthon) multisolitons in the two- 
dimensional (2D) ("pancake-shaped configuration") Bose-Einstein condensate (BEC) with attractive 
interaction. By means of a linear stability analysis, we investigate the stability of these structures 
and show that rotating dipole solitons are stable provided that the number of atoms is small enough. 
The results were confirmed by direct numerical simulations of the 2D Gross-Pitaevskii equation. 
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Localized coherent structures, such as fundamental 
solitons, vortices, nonrotating and rotating niultisoUtons 
are universal objects which appear in many nonlinear 
physical systems and, in particular, in Bose-Einstein 
condensates (EEC's). Stability of these nonlinear struc- 
tures is one of the most important questions because of 
its direct connection with the possibility of experimental 
observation of solitons and vortices. 

Detailed investigations of the stability of localized vor- 
tices in an effectively two-dimensional (2D) trapped BEC 
with a negative scattering length (attractive interaction) 
were performed in Ref. and later extended to the 
three-dimensional case 0] (see also Ref. [4]). While vor- 
tex solitons in attractive BEC are strongly unstable in 
free space, the presence of the trapping potential results 
in existence of stable vortices provided that the number 
of particles does not exceed a threshold value [1, |3) Hi • 

Recently, a novel class of 2D spatially localized vor- 
tices with a spatially modulated phase, the so called az- 
imuthons, was introduced in Ref. Azimuthons rep- 
resent intermediate states between the radially symmet- 
ric vortices and rotating soliton clusters. In contrast to 
the linear vortex phase, the phase of the azimuthon is a 
staircaselike nonlinear function of the polar angle. Vari- 
ous kinds of azimuthons have been shown to be stable in 
media with a nonlocal nonlinear response 0, H, H, Hi [HI ■ 

The aim of this Brief Report is to present nonrotating 
multisoliton (in particular, dipole and quadrupole) and 
rotating multisoliton (azimuthon) structures in the 2D 
BEC with attraction and study their stability by a lin- 
ear stability analysis. We show that, in the presence of 
a confining potential, stable azimuthons also exist for a 
medium with a local cubic attractive nonlinearity. Re- 
sults of the linear stability analysis were confirmed by 
direct numerical simulations of the azimuthon dynamics. 

We consider a condensate which is loaded in an axisym- 
metric with respect to the (x, y) plane harmonic trap, and 
tightly confined in the z direction. The dynamics of the 
condensate is described by the Gross-Pitaevskii equation 



(GPE) 



-|-A + |[f7^(x2+z/) + f7l.^] 



(1) 



where ^'(r,t) is the condensate wave function, a ^ is 
the s-wave scattering length. We assume that the axial 
confinement frequency f2z is much larger than the ra- 
dial one fir (the pancake configuration). Then, the 3D 
equation ([T]) can be reduced to an effective GPE in two 
dimensions [2] (for a detailed discussion of the applicabil- 
ity of the 2D approximation see Ref. 0)- The standard 
reduction procedure [H, leads to the following 2D 
equation 

= -A^^ + ^(x^ + y^)^ - al^l V, (2) 

where appropriate dimensionless units are used, and 
A_L — d/dx^ -\- djdy^, a — ±1, where the +(— ) sign 
corresponds to attractive (repulsive) contact interaction. 

Equation conserves the 2D norm (the normalized 
number of particles) 



dxdy, 



(3) 



and energy 



(4) 

Stationary solutions of Eq. ([2]) in the form 

■ipix, y, t) = Lp{x, y) exp{-ifit), (5) 

where fj, is the chemical potential, resolve the variational 
problem SS — for the functional 



S = E- ^iN. 



(6) 
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Following Ref. f3], let us consider the trial function ip — 
A(j)a(r/a,9), where 4>a is some test function, A and a is 



2 



the amplitude and characteristic width of the stationary 
state, respectively. Then, the functional S takes the form 

S{a, A) = aA^ ~ (3A^a^ + ^-fA^a^ - ii5A^a\ (7) 
where the integral coefficients are 



(9) 



The Euler-Lagrange equations dS/da = and 
dS/dA = for the functional Eq. ([7]) give expressions 
for the width a and amplitude A [5j 



a^,) = + (10) 



^^(m) = 



2/3 



(11) 



The 2D norm N and energy £' then read 



N{ii)=A\p)a\p)5 



^a^(Ai)/r.2„2 



2/3 



(r!^a^7 - 2^^^), (12) 



(13) 



In the following, we consider attractive short-range in- 
teractions and set cr = 1. It then follows from Eq. (fTT|) 
that solutions exist provided that the chemical potential 
^ does not exceed a critical value /x* 



(14) 



To proceed further, we take a test function (j)a in the form 

(t)^{^) = ^'''Li^'\^^)e-^'/^{cosm0 + ipsinme), (15) 

where m is an integer, is the aziniuthal angle, ^ p ^ 1, 
and L^™^ (x) is the nth (i. e. with n zeros) generalized 
Laguerre polynomial. The parameter p determines the 
modulation depth of the soliton intensity. Note that the 
case p = corresponds to the nonrotating multisolitons 
(e. g. m = 1 to a dipole, m = 2 to a quadrupole etc.), 
while the opposite case p = 1 corresponds to the radially 
symmetric vortices. The intermediate case < p < 1 
corresponds to the rotating azimuthons. Since the vor- 
tices and azimuthons have a nontrivial phase, they carry 
out the nonzero z-component of the angular momentum 



= Im J [V'*(r X \'±ip)]^dxdy, (16) 
which can be expressed as 

(17) 



When n 7^ 0, the ansatz Eq. ([T5|) represents a rather 
exotic structure with n nodes. 

Inserting Eq. (fT5|) into Eqs. dH) and ([9]), one can de- 
termine the coefficients a, /3, 7 and 5. For the nodeless 
case (n = 0) we have 

4™^=7r^ = |(p' + l)(- + l)!, (18) 
^(3p4 + 2p2 + 3)(2m)!, (19) 

e(m) TT 



5r = 7;{p' + lW- 



(20) 



For the structures with one node (ri ~ 1), one can obtain 



(„) _ 7r(m + 3) 



(p2 + l)(m+ 1)!, 



(21) 



/3i 



(m) 



7r(3m2 -I- 5m -I- 2) 



5[ 



22m+7 
(m) _ 

~ 2 



(3/ 2/ + 3) (2m)!, (22) 
(p2 + l)(m+l)!, (23) 



7i 



(m) _ 



2(m + 2) 



(p2 + l)(m + 3)!. 



(24) 



Substituting Eqs. (ITHl)-(I2ni) or Eqs. (I1I1)-(E11) into Eqs. 
(fTU)) and pT|) . we get the width and amplitude of the 
corresponding nonlinear structures. In what follows we 
restrict ourselves to the case of the nodeless states {n = 0) 
and, in addition, set JIq = 2. 

The dependences ii{N) and E{N) obtained from Eqs. 
and p3p for different values p are plotted in Fig. [1] 
for m = 1 and m = 2. The curve 4 corresponding radially 
symmetric vortices p — 1 coincides with that obtained in 
Ref. It follows from Eq. ^ that the critical value 

for the chemical potential is /i* = 2(m -I- 1) and does not 
depend on p. The asymptotic values Nmaxi^J■ = —00), 
which determine the (formal) collapse threshold, decrease 
with decreasing p. Results of the variational analysis are 
found to be in good agreement with numerical simula- 
tions (see below). 

Generally speaking, using the relaxation technique 
similar to one described in Ref. (l3 | and choosing an 
appropriate initial guess, one can find numerically vor- 
tex and azimuthon solutions of Eq. ([2]) on Cartesian 
grid. Under this, the parameter p (modulational depth), 
which is similar to the one in Eq. ITSl) , can be introduced 
in the following way: 



p — max |Im ^fj/ max |Re 5'| 



(25) 



However, the choice of initial guess (to achieve conver- 
gence) is extremely difficult and time consuming, and, 
moreover, we were not able to find azimuthon solutions 
with arbitrary p. Instead, we use an approximate but 
much simpler variational approach and introduce the fol- 
lowing ansatz in polar coordinates (r, 9) 0] 



ip{r, t) ~ J7(r)(cos mO -\- ipsinmS 



-ifit 



(26) 



where U (r) is a real function. Inserting the ansatz 
into the action Eq. ([5]), integrating over 6, but keeping 
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FIG. 1: (a) Chemical potential fi and (b) energy E versus 
normalized number of atoms A'' for variational nodeless (n — 
0) solutions Eq. (|15|l with p = (curve 1), p — 0.3 (curve 2), 
p = 0.6 (curve 3), and p = 1 (curve 4) for m — 1; (c) and (d) 
the same for m — 2. 



an arbitrary dependence U{r), one can then obtain the 
corresponding Euler-Lagrange equation, 

(fU IdU { 2 rn^\ rr -.rr^ „ 

TrT + -T^+U^-'^'-^W + /(p)f/'-o, (27) 



dr2 
where 



r dr 



fip) 



3p^ + + 3 
4(p2 + 1) 



(28) 



Equation was solved numerically with boundary con- 
ditions U r'™! at r ^ 0, and U ^ at r ^ oo. In Fig. 
[2] we demonstrate an example of the azimuthon with two 
intensity peaks (i. e. with the topological charge m — 1), 
p = 0.7, and fi = 2.7. For fixed chemical potential /i and 
integer to, there is a family of azimuthon solutions with 
different p. Note that the ansatz (|26p represents only 
particular class of the azimuthons. More general form 
of the azimuthon solutions, which are characterized by 
two independent integer numbers (number of peaks is 
generally independent on the topological charge to) , was 
introduced in Ref. 

To study the stability of the stationary solutions, we 
represent the wave function in the form 



i^{r,t) = [(^o(r)+e(r,t)]e- 



(29) 



where the stationary solution ipo is perturbed by a small 
perturbation e. As usual, one could then take 



(30) 



and consider the corresponding eigenvalue problem, but, 
in our case with p ^ 1, this approach turns out to be 





FIG. 2; (a) Amplitude \ip\ and (b) phase arg ip distributions of 
the azimuthon with two intensity peaks (m = 1) and p — 0.7, 
M = 2.7. 




FIG. 3: The growth rates 7 as functions of the chemical po- 
tential fi for m = 1 and different p. 



ineffective. Indeed, the linearization of Eq. ^ around 
ipo in e leads to the eigenvalue problem 



Tip+ - 



* 2 



ujip+, 



(31) 
(32) 



where T = /i + — {x"^ + y'^) + 2\ipo\'^ and lu are eigenval- 
ues. Nonzero imaginary parts in to imply the instability 
of the state <po(r) with max|Ima;| being the instability 
growth rate. For radially symmetric vortices with p — 1, 
azimuthal perturbations turn out to be the most danger- 
ous, and the problem Eqs. ([51]) and ((5^ can be reduced 
to ID (radial) one and can then be easily solved with 
high accuracy 0. The situation, however, changes dra- 
matically for p ^ 1. In this case the radial symmetry is 
absent, and one must solve the full 2D eigenvalue prob- 
lem. The problem on a, N x N spatial grid implies a 
2A^2 X 2N'^ complex nonsymmetric matrix and, for rea- 
sonable iV (say, N > 200), represents a formidable task. 
Instead, after inserting Eq. P5|) into Eq. we solved 
the Cauchy problem for the linearized equation 



de 



A^e - [x^ + y^)e + 2\ipo\h + ifle* = (33) 



with some initial perturbations e. The final results are 
not sensitive to the specific form of e{x,y,0). If the dy- 
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FIG. 4: Top row - unstable dynamics of the dipole p = with 
fj, = 3.3; bottom row - stable evolution of the azimuthon with 
m = 1, p = 0.7, and p = 2.7. 



namics is unstable, the corresponding solutions e(x, y, t), 
undergoing, generally speaking, oscillations, grow expo- 
nentially in time and an estimate for the growth rate 7 
can be written as 



1 f Pjt + At) 



(34) 



where P{t) — J \e\'^dxdy, t and At are assumed to be 
large enough. 

In Fig. [3] we plot the growth rates 7 as functions of 
the chemical potential fi for m = 1 and several different 
values of p. The case p = 1 corresponds to radially sym- 



metric vortices with the topological charge m = I and 
was investigated in detail in Ref. Q. The corresponding 
curve in Fig. [3] coincides with that obtained in Ref. 0] ■ 
The linear stability analysis shows that for solutions with 
m — 1 and p < 0.7, the growth rate of perturbations 
7 7^ for all fi so that there is no the stability region. 
The situation, however, changes when p > 0.7. Under 
this, the growth rate of perturbations 7 = if the chem- 
ical potential exceeds some critical value /ic. The sta- 
bility window appears and the azimuthons with p > 0.7 
and fXc < IJ- < The critical value /ic monotonically 
increases from ~ 2.45 (for p = 0.7) to ^ 2.6 (for 
p = 1). All solutions with m — 2 turn out to be unstable 
(for p = 1 it was shown in Ref. 0). 

To verify the results of the linear analysis, we solved 
numerically the dynamical equation ([2]) initialized with 
our computed solutions with added Gaussian noise. The 
initial condition was taken in the form ipo[l + z^^(a;,y)], 
where ipQ{x,y) is the numerically calculated solution, 
^{x,y) is the white gaussian noise with variance = 1 
and the parameter of perturbation 1/ = 0.005 0.1. The 
unstable dynamics of the nonrotating dipole (p = 0) with 
fi = 3.3 is illustrated in the top row of FigHl Stable evo- 
lution of the azimuthon with p — 0.7 and /i = 2.7 (i. e. 
in the region of the stability) is shown in the bottom row 
of FigHl The azimuthon cleans up itself from the noise 
and survives over hundreds of the rotational periods. 

In conclusion, we have presented nonrotating and ro- 
tating (azimuthon) multisolitons in an effectively 2D 
(" pancake-shaped configuration" ) Bose- Einstein conden- 
sate with attractive interaction and parabolic trapping 
potential. We have performed a linear stability analysis 
of these structures and demonstrated that azimuthons 
with two intensity peaks (rotating dipoles) and with not 
too small modulational depth can be stable if the number 
of particles is below some critical value. The nonrotat- 
ing multisolitons (dipoles and all high-order multipoles) 
appear to be unstable. 
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